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Abstract 

A Bayesian analysis of the probability of a signal in the presence 
of background is developed, and criteria are proposed for claiming ev- 
idence for, or the discovery of a signal. The method is general and in 
particular applicable to sparsely populated spectra. Monte Carlo tech- 
niques to evaluate the sensitivity of an experiment are described. 

As an example, the method is used to calculate the sensitivity of the 
GERDA experiment to neutrinoless double beta decay. 

1 Introduction 

In the analysis of sparsely populated spectra common approximations, valid 
only for large numbers, fail for the small number of events encountered. A 
Bayesian analysis of the probability of a signal in the presence of background is 
developed, and criteria are proposed for claiming evidence for, or the discovery 
of a signal. It is independent of the physics case and can be applied to a variety 
of situations. 

To make predictions about possible outcomes of an experiment, distribu- 
tions of quantities under study are calculated. As an approximation, ensem- 
bles, sets of Monte Carlo data which mimic the expected spectrum, are ran- 
domly generated and analyzed. The frequency distributions of output param- 
eters of the Bayesian analysis are interpreted as probability densities and are 
used to evaluate the sensitivity of the experiment to the process under study. 

As an example, the analysis method is used to estimate the sensitivity of 
the GERDA experiment P to neutrinoless double beta decay. 

The analysis strategy is introduced in section [21 The generation of ensem- 
bles and the application of the method onto those is discussed in section El 
The application of the analysis method in the GERDA experiment is given as 
an example in section 01 where the sensitivity of the experiment is evaluated. 



2 Spectral analysis 

A common situation in the analysis of data is the following: two types of 
processes (referred to as signal and background in the following) potentially 
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contribute to a measured spectrum. The basic questions which are to be 
answered can be phrased as: What is the contribution of the signal process 
to the observed spectrum? What is the probability that the spectrum is due to 
background only? Given a model for the signal and background, what is the 
(most probable) parameter value describing the number of signal events in the 
spectrum? In case no signal is observed, what is the limit that can be set on the 
signal contribution? The analysis method introduced in this paper is based on 
Bayes' Theorem and developed to answer these questions and is in particular 
suitable for spectra with a small number of events. 
The assumptions for the analysis are 

• The spectrum is confined to a certain region of interest. 

• The spectral shape of a possible signal is known. 

• The spectral shape of the background is known^. 

• The spectrum is divided into bins and the event numbers in the bins 
follow Poisson distributions. 

The analysis consists of two steps. First, the probability that the observed 
spectrum is due to background only is calculated. If this probability is less 
then an a priori defined value, the discovery (or evidence) criterion, the signal 
process is assumed to contribute to the spectrum and a discovery (or evidence) 
is claimed. If the process is known to exist, this step is skipped. Based on 
the outcome, in a second step the signal contribution is either estimated or an 
upper limit for the signal contribution is calculated. 

2.1 Hypothesis test 

In the following, H denotes the hypothesis that the observed spectrum is due 
to background only; the negation, interpreted here as the hypothesis that 
the signal process contributes to the spectrum^, is labeled H. The conditional 
probabilities for the hypotheses H and H to be true or not, given the measured 
spectrum are labeled p(if|spectrum) and p(if| spectrum), respectively. They 
obey the following relation: 

p(if I spectrum) + p(i7| spectrum) = 1 . (1) 



The conditional probabilities for H and H can be calculated using Bayes' 
Theorem |2|: 



This assumption and the previous can be removed in a straightforward way with the 
introduction of additional prior densities. 

^ Since the shape of the background spectrum is assumed to be known the case of unknown 
background sources contributing to the measured spectrum is ignored. However, the overall 
level of background is allowed to vary. 




p(i!/" I spectrum) 



^(spectrumlif) ■ Pq{H\I) 



(2) 



p(spectrum) 
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and 



p(if|spectrum) 



p(spectrum|i7) • po{H\I) 
p(spectrum) 



(3) 



where p(spectrum|if) and p (spectrum | if) are the conditional probabihties to 
find the observed spectrum given that the hypothesis H is true or not true, 
respectively and po{H\I) and po{H\I) are the prior probabilities for H and H. 
The values of po{H\I) and po{H\I) are chosen depending on additional infor- 
mation, /, such as existing knowledge from previous experiments and model 
predictions. In the following, the symbol / is dropped but it should be under- 
stood that all available information is used in the evaluation of probabilities. 
The probability p(spectrum) is rewritten as 

p(spectrum) = p (spectrum ■ po{H) -|- p (spectrum | if) ■ po{H) (4) 

The probabilities p(spectrum|ii) and p(spectrum|i7) can be decomposed 
in terms of the expected number of signal events, S, and the expected number 
of background events, B: 

p(spectrum|ii) = J p{spectrum\B) ■ pq{B) dB, (5) 

p(spectrum|i7) = J p{spectrum\S, B) ■ pq{S) • po{B) dS dB , (6) 

where p(spectrum|i?) and p(spcctrum|S', B) are the conditional probabilities 
to obtain the measured spectrum. Further, po{S) and po{B) are the prior prob- 
abilities for the number of signal and background events, respectively. They 
are assumed to be uncorrelated, and are chosen depending on the knowledge 
from previous experiments, supporting measurements and models. 

The observed number of events in the ith bin of the spectrum is denoted 
rii. Assuming the fluctuations in the bins of the spectrum to be uncorrelated 
the probability to observe the measured spectrum, given B (in case H is true) 
or the set S, B (in case H is true), is simply the product of the probabilities 
to observe the values, {rii}. The expected number of events in the ith bin, 
Aj, can be expressed in terms of S and B: 



where fs{E) and /b(-E') are the normahzed shapes of the known signal and 

background spectra, respectively, and AEi is the width of the zth bin. The 
letter E suggests an energy bin, but the binning can be performed in any 
quantity of interest. The number of events in each bin can fluctuate around 



Aj = K{S, B) 

^ S - [ fs{E) dE + B- [ h{E) dE, 



(7) 
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Aj according to a Poisson distribution. This yields 



N 



p(spectrum|i?) 



i=l 

N 



p(spectrum|5, S) = ^r^ -MS, b) 



1=1 



nA 



(8) 



(9) 



In summary, the probabihty for H to be true, given the measured spectrum. 



is: 



I spectrum) 



(10) 





/n|re-"'-Po(-B) dB 




r \"i 1 
!W^^-^-Po{B) dB 




r \"i 1 
!Y{^e-^-Po{B)vo{S) dB dS 


■Po{H) 



with Aj calculated according to ((Tj). Evidence for a signal or a discovery can be 
decided based on the resulting value for p(if| spectrum). It should be empha- 
sized that the discovery criterion has to be chosen before the data is analyzed. 
A value of p(if|spectrum) < 0.0001 is proposed for the discovery criterion, 
whereas a value of p(if|spectrum) < 0.01 can be considered to give evidence 
for H. 

The analysis can be easily extended to include uncertainties in the knowl- 
edge of relevant quantities. For example, if the spectrum is plotted as a func- 
tion of energy, and the energy scale has an uncertainty, then equations ()8p9|l 
can be rewritten as 



p(spectrum|i?) 
p(spectrum|5', B) 



N 



j-j A,(0, i^|fc)"' ^_,,^(o, Bik) 



i=l 
■ N 



i=l 



Po{k)dk (11) 
Po{k)dk . (12) 



where Xi{S, B\k) is the expected number of events for a given energy scale 
factor k andpo(^) is the probability density for k (e.g., a Gaussian distribution 
centered on = 1). 



2.2 Signal parameter estimate 

In case the spectrum fulfills the requirement of evidence or discovery, the 
number of signal events can be estimated from the data. The probability 
that the observed spectrum can be explained by the set of parameters S and 
B, making again use of Bayes' Theorem, is: 



p{S, -B I spectrum) 



p(spectrum|S', B) ■ pq{S) ■ po{B) 
f p(spectrum|S', B) ■ Po{S) ■ po{B) dS dB 



(13) 
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In order to estimate the signal contribution the probabihty p{S, 5 1 spectrum) 
is marginahzed with respect to B: 

p(5'|spectrum) = J p(S, -B|spectrum) dB . (14) 

The mode of this distribution, S*, i.e., the value of S which maximizes 
p( 5*1 spectrum), can be used as an estimator for the signal contribution. The 
standard uncertainty on S can be evaluated from 

rSi6 

/ p( 5*1 spectrum) dS = 0.16 
Jo 

/ p(S' I spectrum) dS = 0.84 
Jo 



such that the results can be quoted as 

-(584-^*) 



2.3 Setting limits on the signal parameter 

In case the requirement for an observation of the signal process is not fulfilled 
an upper limit on the number of signal events is calculated. For example, a 
90% probability lower limit is calculated by integrating Equation (|T^ to 90% 
probability: 

/"S'go 

/ p(^lspectrum) dS = 0.90 . (16) 
Jo 

5*90 is the 90% probability upper limit on the number of signal events. It 
should be noted that in this case it is assumed that H is true but the signal 
process is too weak to significantly contribute to the spectrum. 



3 Making predictions - ensemble tests 

In order to predict the outcome of an experiment distributions of the quantities 
under study can be calculated. This is done numerically by generating possible 
spectra and subsequently analyzing these. The spectra are typically generated 
from Monte Carlo simulations of signal and background events. For a given 
ensemble, the number of signal and background events, 5*0 and Bq are fixed 
and a random number of events are collected according to Poisson distributions 
with means 5*0 and Bq. From each ensemble a spectrum is extracted and the 
analysis described above is applied. The analysis chain is shown in Figure ^ 
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The output parameters, such as the conditional probabihty for H, 
I spectrum), are histogrammed and the frequency distribution is interpreted 
as the probabihty density for the parameter under study. As examples, the 
mean value and the 16% to 84% probability intervals can be deduced and used 
to predict the outcome of the experiment. This approach is referred to as 
ensemble tests. 

Systematic uncertainties, such as the influence of energy resolution, miscal- 
ibration or signal and background efficiencies, can be estimated by analyzing 
ensembles which are generated under different assumptions. 




_ ^ Ensemble 1 } Speclrum 1 } '- (p(H|spectrum 1) ) 

/ [Ensemble~2} - * (Speclrum 2 } "- (plHlspectrum 2) ) 

(Ensemble's) - " (Spectrum 3 } '- (p(H|spectrum 3) ) 

Analysis^ 
Sets of spectra- 



Ensembles: 




sets of events 






[Ensemble 4} - {Spectrum 4 ) - - (p(H|speclrum 4) ) 

[Ensemble ..} - - {Spectrum ...) - [p(H|spectrum ...)) 

[EnsembleTi} - - {Spectrum n } - - [p(H|spectrum n) ] 



Figure 1: Analysis chain. The Monte Carlo generator (MC) generates a pool 
which consists of signal and background events. An ensemble is defined as a 
set of events representing a possible outcome of an experiment. The number 
of events are randomly chosen according to the parameters 5*0 and Bq. From 
each ensemble a spectrum is extracted and subsequently analyzed. The prob- 
ability j9(iir I spectrum) for each spectrum is depicted here as the outcome of 
the analysis. 



4 Sensitivity of the GERDA experiment 

In the following, the GERDA experiment is introduced and the Bayesian 
analysis method, developed in section |2l is applied on Monte Carlo data in 
order to predict possible outcomes of the experiment. 

4.1 Neutrinoless double beta decay and the GERDA ex- 
periment 

The GERmanium Detector Array, GERDA |T], is a new experiment to search 
for neutrinoless double beta decay (0z//3/9) of the germanium isotope ^^Ge. 
Neutrinoless double beta decay is a second order weak process which is pre- 
dicted to occur if the neutrino is a Majorana particle. The half-life of the 
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process is a function of the neutrino masses, their mixing angles, and the CP 
phases. Today, 90% C.L. hmits on the half-hfe for neutrinoless double beta 
decay of ^^Ge exist and come from the Heidelberg-Moscow IH] and IGEX |3] 
experiments. They are T1/2 > 1-9 ■ 10^^ years and T1/2 > 1.6 ■ 10^^ years, 
respectively. A positive claim was given by parts of the Heidelberg-Moscow 
collaboration with a 3 a range of T1/2 = (0.7 — 4.2) ■ 10^^ years and a best value 
of Ti/2 = 1.2 ■ 10^5 years 0. 

A total exposure (measured in kg-years of operating the germanium diodes) 
of at least 100 kg-years should be collected during the run-time of the GERDA 
experiment. The germanium diodes are enriched in the isotope ^^Ge to a level 
of about 86%. One of the most ambitious goals of the experiment is the 
envisioned background level of 10~^ counts/ (kg-keV-y). This is two orders of 
magnitude below the background index observed in previous experiments |HIIE]- 
For an exposure of 100 kg-years the expected number of background events in 
the 100 keV wide region of interest is approximately 10. Using the present 
best limit on the half-life less than 20 Oz//5/5-events are expected within a much 
smaller window. The number of expected Oi//5/5-events, Sq, is correlated with 
the half-life of the process via 

5o^ln2-«:-M-6sig-^--^, (17) 

Ma J 1/2 

where k = 0.86 is the enrichment factor, M is the mass of germanium in grams, 
TVa is Avogadro's constant and t is the measuring time. Ma is the atomic mass 
and esig is the signal efficiency, estimated from Monte Carlo data to be 87%. 



4.2 Expected spectral shapes and prior probabilities 

In GERDA, the energy spectrum in the region around 2 MeV is expected 
to be populated by events from various background processes. The signature 
of neutrinoless double beta decay, the signal process, is a sharp spectral line 
at the Q^^- value which for the germanium isotope ^'^Ge is 2 039 keV. In the 
following, the region of interest is defined as an energy window of ±50 keV 
around the Q/3/3-value. The shape of the background spectrum is assumed to 
be fiat, i.e. fB{E) = const. The shape of the signal contribution is assumed to 
be Gaussian with a mean value at the Q/3/3-value. The energy resolution of the 
germanium detectors in the GERDA setup is expected to be 5 keV (FWHM), 
corresponding to a width of the signal Gaussian of a ~ 2.1 keV. 

For the calculation of the sensitivity, ensembles are generated according 
to (1) the exposure, (2) the half-life of the Oi//3/3-process which is translated 
into the number of expected signal events, Sq, in the spectrum, and (2) the 
background index in the region of interest which is translated into the number 
of expected background events, Bq. The number of signal and background 
events in each ensemble fiuctuate around their expectation values 5*0 and Bq 
according to a Poisson distribution. For each set of input parameters 1000 
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ensembles are generated. An energy spectrum is extracted from each ensemble 
with a bin size of 1 keV. 

In order to calculate the probability that the spectrum is due to background 
processes only, the prior probabilities for the hypothesis H and H have to be 
fixed, as well as those for the signal and background contributions. This is 
a key step in the Bayesian analysis. Given the lack of theoretical consensus 
on the Majorana nature of neutrinos and the cloudy experimental picture, the 
prior probabilities for H and H are chosen to be equal, i.e. 

p,{H) = 0.5, (18) 
Po(H) = 0.5. (19) 

The prior probability for the number of expected signal events, assuming 
H, is assumed fiat up to a maximum value, Smax, consistent with existing 
limits^. It should be noted that the setting of the prior probability for H is 
dependent on the maximum allowed signal rate. S^ax was chosen in such a 
way that the probability for the hypothesis H is reasonably assumed to be 
50 %. The effect of choosing a different prior for the number of signal events 
is discussed below. 

The background contribution B is assumed to be known within some un- 
certainty (recall that the shape of the background is however fixed). The prior 
probability for B is chosen to be Gaussian with mean value /xb = -Bo aiid width 
(Jb = Bq/2. The prior probabilities for the expected signal and background 
contributions are 

Po{S) = -J—, < S < Smax, Po{S) = otherwise, (20) 

*^max 

(-B-mb)2 

PoiB) = 'I"^^,2 , B>0, poiB) = 0B<0. (21) 



Jo 



4.3 Examples 

As an example. Figure |2l (top, left) shows a spectrum from Monte Carlo data 
generated under the assumptions of a half-life of 2 ■ 10^^ years, a background 
index of 10~^ counts/(kg-keV-y) and an exposure of 100 kg-years. This cor- 
responds to 5*0 = 20.5 and Bq = 10.0. The (20) signal and (8) background 
events are indicated by a solid and dashed line, respectively. Figure 12 (top, 
right) shows p(S'|spectrum) for the same spectrum. The mode of the distribu- 
tion is S* = 19.8, consistent with the number of signal events in the spectrum. 
Figure El (bottom, left) shows the distribution of S* for 1000 ensembles gen- 
erated under the same assumptions. The average number of S* = 20.3, in 
agreement with the average number of generated signal events, 20.4. Figure El 
(bottom, right) shows the distribution of the logj9(-ff|spectrum) for ensembles 

^Sniax was calculated using Eg uat ion 1 1 71 assuming a half-life of T'1/2 = 0.5 • 10^^ years. 
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generated under the same assumptions. More than 97% of the ensembles have 
a probabihty spectrum) of less than 0.01%. I.e., a discovery could not be 
claimed for less than 3% of experiments under these conditions. 




5 10 15 20 25 30 35 40 -40 -35 -30 -25 -20 -15 -10 -5 

S* log p(H|spectrum, I) 



Figure 2: The spectrum (top, left) was randomly generated under the assump- 
tions of a half-life of 2-10^^ years, a background index of 10~^ counts/ (kg-keV-y) 
and an exposure of 100 kg-years. The signal events are indicated by a solid 
line, the background events by a dashed line. The probability density for S 
(top, right) peaks at 19.8 which is consistent with the number of signal events, 
20, in the spectrum. The distribution of the estimated number of signal events 
(bottom, left) as well as the distribution of the logp(if|spectrum) (bottom, 
right) are calculated from ensembles generated under the same assumptions. 

In order to simulate the case in which only lower limits on the half-life 
of the Oi//3/5-process are set, ensembles are generated without signal contribu- 
tion, i.e. 5*0 = 0. As an example. Fig. IHl (top, left) shows a spectrum from 
Monte Carlo data generated under the assumptions of a background index of 
10"'^ counts/ (kg-keV-y) and an exposure of 100 kg-years. No signal events are 
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present in the spectrum. 

Figure 01 (top, right) shows the marginahzed probabihty density for S, 
p(S'lspectrum), for the same spectrum. The mode of is events. 

Figure El (bottom, left) shows the distribution of the hmit (90% probabil- 
ity) of the signal contribution for 1000 ensembles generated under the same 
assumptions. The average limit is 3.1. 

Figure El (bottom, right) shows the distribution of the p(if|spectrum) for 
ensembles generated under the same assumptions. For none of the ensembles 
could a discovery be claimed. 




2000 2020 2040 2060 2080 5 10 15 20 



E [keV] s 




$800 



2 4 6 8 10 

S (90% C.L. upper limit) 




0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 

p(l-l|spectrum, I) 



Figure 3: The spectrum (top, left) was randomly generated under the assump- 
tions of a background index of 10^^ counts/ (kg-keV-y) and an exposure of 
100 kg-years. No signal events are present in the spectrum. The probability 
density for S for the same spectrum (top, right) peaks at 0. The distribution 
of the limit (90% probability) of the signal contribution (bottom, left) as well 
as the distribution of the p(if| spectrum) (bottom, right) are calculated from 
ensembles generated under the same assumptions. 
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4.4 Sensitivity 

For the ensembles generated without signal contribution the mean of the 90% 
probability lower limit on the half-life is shown in Figure lU as a function 
of the exposure for different background indices. In case no background is 
present the limit scales linearly with the exposure. With increasing back- 
ground contribution the limit on the half-life increases more slowly. For the 
envisioned background index of 10"^ counts/(kg-keV-y) and an expected ex- 
posure of 100 kg-years an average lower limit of T1/2 > 13.5 ■ 10^^ years can be 
set. For the same exposure, the average lower limit is T1/2 > 6.0 ■ 10^^ years 
and Ti/2 > 18.5 ■ 10^^ years for background indices of 10~^ counts/ (kg-keV-y) 
and 10"'^ counts/ (kg-keV-y), respectively. 

Using the nuclear matrix elements quoted in ^ the lower limit on the half- 
life of the Ou (3 (3-pTocess can be translated into an upper limit on the effective 
Majorana neutrino mass, {mpp), via 

{m,,) = ■ G^n-'^' ■ (22) 

where G^'^ is a phase space factor and {M^'^) is the nuclear matrix element. 
Figure 13 also shows the expected 90% probability upper limit on the effective 
Majorana neutrino mass as a function of the exposure. With a background 
index of 10~^ counts/(kg-keV-y) and an exposure of 100 kg-years, an upper 
limit of {rrtpp) < 200 meV could be set assuming no Oi//3/5-events are observed. 

Figure El shows the half-life for which 50% of the experiments would re- 
port a discovery of neutrinoless double beta decay as a function of the expo- 
sure for different background indices. For the envisioned background index of 
10~^ counts/ (kg-keV-y) and an expected exposure of 100 kg-years this half-life 
is 5 - 10^^ years. 

Using the same matrix elements from reference [Zj, the half-life is trans- 
formed into an effective Majorana neutrino mass. The mass for which 50% 
of the experiments would report a discovery is shown in Figure El (bottom) 
as a function of the exposure and for different background conditions. For an 
exposure of 100 kg-years and a background index of 10~^ counts/ (kg-keV-y) 
neutrinoless double beta decay could be discovered for an effective Majorana 
neutrino mass of 350 meV (with a 50% probability). 
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50 100 150 200 
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Figure 4: The upper plot shows the expected 90% probabihty lower limit 
on the half-life for neutrinoless double beta decay versus the exposure under 
different background conditions. Also shown is the half-life for the claimed 
observation The lower plot shows the expected 90% probability upper limit 
on the effective Majorana neutrino mass versus the exposure under different 
background conditions. The effective Majorana neutrino mass for the claimed 
observation is also shown. All mass values were determined from the half-life 
using the matrix element reported in 
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Figure 5: Top: the half-life for which 50% of the experiments would report a 
discovery, i.e. a probability that the spectrum is due to background processes 
only, p(if|spectrum), of less than 0.01%, is plotted versus the exposure under 
different background conditions. Bottom: the effective Majorana neutrino 
mass for which 50% of the experiments would report a discovery versus the 
exposure under different background conditions. The mass was determined 
from the half-life using the matrix element reported in [7j. 
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4.5 Influence of the prior probabilities 



In order to study the influence of tlie prior probabilities on the outcome of 
the experiment, the prior probability for the number of expected signal events, 
Pq{S), was varied. Three different prior probabilities were studied: 

• flat prior: po{S) oc const. , 

• pessimistic prior: po{S) oc e~^^^^ , 

• peaking prior: po{S) oc e^~^/^ /S'^, 

where S is the number of events corresponding to a half-life of 1.2-10^^ years 
and S < Smax- For a background index of 10~^ counts/(kg keV y) and an ex- 
posure of 100 kg years the limit strongly depends on the chosen prior. For the 
pessimistic prior probability the limit which can be set on the half-life is about 
10% higher than that for the flat prior probability. In comparison, the peaking 
prior gives a 50% lower limit compared to the flat prior. This study makes 
the role of priors clear. If an opinion is initially strongly held, then substantial 
data is needed to change it. In the scientific context, consensus priors should 
be strived for. 



5 Conclusions 

An analysis method, based on Bayes' Theorem, was developed which can be 
used to evaluate the probability that a spectrum can be explained by back- 
ground processes alone, and thereby determine whether a signal process is 
necessary. A criterion for claiming evidence for, or discovery of, a signal was 
proposed. Monte Carlo techniques were described to make predictions about 
the possible outcomes of the experiments and to evaluate the sensitivity for 
the process under study. 

As an example the method was applied to the case of the GERDA ex- 
periment for which the sensitivity to neutrinoless double beta decay of ^^Ge 
was calculated. With a background index of 10~^ counts/ (kg-keV-y) and an 
exposure of 100 kg-years the sensitivity of the half-hfe of the Oi//3/3-process is 
expected to be 13.5 • 10^^ years. 
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